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ABSTRACT 

The non-orthogonality of algebraic polynomials of field coordinates tradition- 
ally used to model field-dependent corrections to astrometric measurements, gives 
rise to subtle adverse effects. In particular, certain field dependent perturbations 
in the observational data propagate into the adjusted coefficients with consid- 
erable magnification. We explain how the worst perturbation, resulting in the 
largest solution error, can be computed for a given non-orthogonal distortion 
model. An algebraic distortion model of full rank can be converted into a fully 
orthonormal model based on the Zernike polynomials for a circular field of view, 
or a basis of functions constructed from the original model by a variant of the 
Gram-Schmidt orthogonalization process for a rectangular field of view. The 
relative significance of orthonormal distortion terms is assessed simply by the 
numerical values of the corresponding coefficients. Orthonormal distortion mod- 
els are easily extendable when the distribution of residuals indicate the presence 
of higher order terms. 

Subject headings: astrometry — methods: data analysis 

1. Introduction 

The need to exactly characterize the mapping of celestial angular coordinates onto 
Cartesian coordinates in the plane of a detector, expressed in suitable units, arises in many 
astronomical observations. Often called distortion models, such mapping models usually 
include the basic geometry of an optical projection, possible deviations in the focal length, 
optical aberrations, etc. The propagation of accidental error in traditional distortion models 
based on algebraic polynomials of field coordinates has been studied in depth, and the 
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use of orthogonal functions has been advocated in the astrometric hterature. However, it 
appears that simple algebraic polynom ials are still in us e for la rge- volume and high-precision 
data reduction systems. For example, iMcArthur et al.l (|2005[ ) represent the field-dependent 
errors in the FGS observations with the HST with two independent fifth order polynomials 
for X and a total of 16 terms in each polynomial. The distortion model is expressed 
through 32 coefficients, which were initially determined by ray-traces on the ground, but 
later updated several times by means of special calibration observations of M35 while in- 
orbit. The high- volume astrometric data reductions in the Sloan Digital Sky Survey (SDSS, 
Pier et al.l 120031 ) were performed with a two-step characterization of field distortions. For 
each CCD separately, the geometric part of the coordinate transformation was modeled 
as a third order polynomial of y only, because the observations were taken in the time- 
delayed integrati on mode. This t ransfo rmation was further superimposed with a 2D affine 
transformation. iBellini fc BedinI ( 120091 ) used independent 3rd order algebraic polynomials 



for Ax and Ay to fit field- dependent corrections to a nominal distortion model (sometimes 
called astrometric transfer function), represented by a conformal transformation of field 
coordinates. 

We investigate the propagation of systematic error and conclude that the use of non- 
orthogonal terms can lead to significant errors in the instrument characterization. The 
worst systematic perturbation in the data, producing the largest possible error in the model 
parameters, can be exactly computed for any particular mo del, which is exemplifi ed by the 
JMAPS astrometric project currently under development ( iHennessy et al.ll2010l ). For the 
original JMAPS focal plane model, the worst perturbation is amplified by a factor of 7 with 
respect to the more benign perturbations of the same norm. To avoid this potentially harmful 
build-up of systematic error, we suggest using orthogonal functions on the unit square (for 
a square detector) and explain how an existing non-orthogonal model can be orthogonalized 
through the Gram-Schmidt process. Determining the significance and temporal stability 
of individual model terms, which are nearly independent on sufficiently dense star fields, 
becomes straightforward. 



2. Polynomial plate models 

Traditionally, the transformation of standard coordinates ^, r] (related to the local angu- 
lar coordinates on the celestial sphere through a simple gnomonic projection) into Cartesian 
coordinates in the detector plane x, y is represented by a simple algebraic polynomial in 
powers of x and y. The repercussions of this choice for the propagation of accidental mea- 
surement errors into the astrometric plate reductions have been discussed in numerous papers 
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[e.g. lEichhorrull957l : |Jefferyslll963l : lEichhorn fc Williamall963l : Ide Vegtlll967l ). The main mo- 
tivation for selecting this functional form was probably its mathematical simplicity, but also 
the fact that some common kinds of distortions, such as rotation, scale and tilt are repre- 
sented by monomials in field coordinates. In the following, we will use a model currently 
considered for the JMAPS astrometric satellite, but the techniques developed in this paper 
can be applied to any variety of polynomial distortion models. The JMAPS model is: 



^ = gqx + aiu + 02 + OeX^ + a-jxy + Ogr^ + oioxr^ (1) 
■q = a^x + a^y + as + a^xy + a-jy'^ + agr^ + aioyr'^ 

with r denoting the radius, r = a/xM-^, for simplicity. This model includes the shifts in x 
and y (model parameters a2 and as), rotation, shear and scale (ao through as), tilt (ae and 
ay), cubic radial distortion (aio) and differential radial quadratic distortion ( as and ag). It 
is convenient to represent this 2D model in vector and matrix forms, by stacking all the x 
and y measurements in a 2 x n column vector (where n is the number of reference stars), 
and likewise the standard coordinates ^ and 1], and arranging the model term values in a 2n 
by 11 matrix, so that the least-squares problem takes the form 

" X y 1 xy xr^ 

xy y^ yr^ x y 1 

where 

a = [ ao ai a2 ag a? ag aio a^ a^ as ag ] . (3) 

The system of linear equations on a discrete set of data points {C,i,rii) is solved by the least 
squares method yielding a solution 

a = (F^ F)-^ F^ 

We are concerned with the propagation of the measurement error in the least squares solution. 
Certainly, it depends on the distribution of data points in the focal plane, among other 
factors. Reference stars in the focal plane can be non-uniformly distributed, for example, 
clustered in the center of the field, or in one of the corners, giving rise to large errors in 
some of the parameters, and to enhanced covariances between distortion coefficients. In 
space astrometry projects, however, a single set of distortion coefficients is determined for a 
large collection of frames, representing a large sample of random stellar configurations and 
accruing a sufficiently high number density of reference stars. It is therefore justified for 
practical and theoretical reasons to consider an idealized sampling of the focal plane with a 
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large and dense regular grid of data points. The inner product of two vectors in that case is 
replaced with its form for a Hilbert function space. For example, the inner product of the 
fourth and the sixth columns of matrix F in Eq. [2] is defined as 



(f4,f6 



dx 



56 
45' 



(5) 



We assumed in this equation (without a loss of generality) that the detector area is square 
and the field coordiantes are normalized to unity. 

In practice, when an astrometric instrument is calibrated by its own observations, it is 
convenient to simplify the transformations between the angular coordinates (^, rf) and the 
detector coordinates (x, y) to a simple form, e.g., a gnomonic projection plus a nominal scale 
conversion. The corrections to the nominal parameters and the higher order terms, which 
may be time-dependent, are sought for as small corrections to the detector coordinates, 
{Ax, Ay). The distortion coefficients are small enough to use the same model (such as in 
Eq. [T]) in it s differential form, by repl acing and rj on the left-hand side with Ax and Ay, 
respectively ( lHiltnerlll962t lGreenlll985[ ). Once the small corrections to the nominal distortion 
model are known, both forward and backward transformations between the field coordinates 
can be easily performed. The JMAPS data reduction system is using this approach too, but 
we retain the notations of Eq. [T] in this paper for simplicity. 



3. Error propagation 

An error in the right-hand part of Eq. [2] (i.e., measurement error e) propagates linearly 
into the least-squares solution a, so that 

Aa = F'^e (6) 

where F''' is the pseudoinverse of F. Therefore, the error in the model parameters is defined 
not only by the measurement error e, but also by the properties of the design matrix, i.e., the 
composition of distortion terms. Modeling field-dependent distortio ns becomes an impor- 



tant i ssue defining to some degree the resulting astrometric accuracy. lEichhorn fc Williams 



(119631 ) derived the covariances of plate constants for a series of commonly used models for 
the case of random uncorrelated errors e. They realized that the uncertainty of field charac- 
terization in the presence of random error may be a complex, but quite predictable, function 
of field coordinates. The uncertainties tend to become larger and more complex for more 
sophisticated polynomial models with higher order terms. They also showed that the resid- 
ual field-dependent error (which they called systematic error) can be larger than the gain in 
precision for unnecessarily involved distortion models. 
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Correlations between some of the terms of a polynomial model are the main reason for 
this potentially harmful feature. A large correlation between two distortion terms should 
always be a concern, because the commonly used performance statistics based on post-fit 
residuals, such as reduced x^, do not capture the intrinsic statistical dependencies within the 
solution. Astrometrists rarely know their instruments well enough in advance, to the extent 
that a priori distortion models can be used without some experimentation or verification 
on real data. Even then, it is not clear where to stop in building up a distortion model, 
because each additional term usually results in smaller residuals, and the solution value can 
be deceptively large, indicating a significant term, being in fact correlated with a number of 
other terms. 

For these reasons, the use of orth ogonal functioii s inste ad of algebrai c polynoni i als ha s 
long been advocated in the literature. iBrosche et al.l (Il989[ ). followed by iBienaymg (119931 ) . 
explained the adverse effects of correlated terms, and considered the ease of determination 
of a minimal-sufficient model to be the main advantage of orthogonalized designs. Having 
determined the level of statistically significant reduction in x^, one can, for example, remove 
orthogonal terms from the model in turns, until this level is achieved with the smallest 
number of terms. 

The effects of systematic errors have largely been unheeded in the literature. By a 
systematic error we understand any deterministic perturbation of the data which can be 
expressed as a function of field coordinates. Obviously, there is an infinite set of possible 
systematic errors, and analysis of their propagation appears to be intractable. Beyond the 
common systematic errors, such as a focal length variation, these effects are hard to pre- 
dict in complex astrometric instruments. However, instead of guessing possible systematic 
errors and simulating their effects, one can pose this problem: is there a specific systematic 
error which propagates into the least squares solution with the largest gain? An astonish- 
i ngly simple answer is obt ained through the Singular Value Decomposition (SVD) technique 
( iGolub fc van LoanI 119961 ). Let 

F = U S (7) 
be the SVD of the design matrix F. The LS solution is 



a = V 



V 



For an overdetermined system with k model terms (i.e., of rank k), is a diagonal matrix 
with k nonzero (positive) numbers in the diagonal in increasing order. Therefore, in the 
projection of the data vector onto the basis, U""" rj]^, only the first k columns of the 
orthogonal matrix U contribute to the solution. In other words, the first k columns of U 
define the basis of the subspace spanned by F and the remainder is the basis of its null 
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spacqll. The LS solution is impervious to any perturbation of data in the null space. If e is 
a systematic perturbation of the data, the squared norm of the resulting systematic error of 
the solution is 

\\5r = 6^6 = el^-^eu (9) 

where eu = U^e is the projection of e onto the basis of F. Recalling that S = diag([cri (72 . . . crkf^) 
and o"i > o"2 > • • • > o'k, the maximum is achieved when e is aligned with = U{:, k). 
This follows from the constrained Lagrange problem, where the Lagrange multiplier is one 
of the (Tj~^ and among all e of unit norm the largest solution error happens when eu = e, 
hence, e = Uk. 
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Fig. 1. — The worst field distortion for the nominal JMAPS distortion, which yields the 
largest error of the focal plane parameters. 



^In Matlab, the significant part of SVD is realized through the economic SVD: [u, sigma, v] =svd(f ,0) . 
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Thus, the kth long singular vector, which is the least signiGcant principal component of 
model F, is the worst systematic error among all possible systematic perturbations of unit 
norm. 



4. The case of JMAPS 

Using a dense grid of points on the unit square, the SVD of the nominal JMAPS 
distortion model (Eq. [2]) can be computed and the worst systematic error numerically 
derived. The result is a higher order field distortion shown in Fig. [TJ It is similar to the 
scale in the center of the field of view, but instead of growing linearly with radius, it turns 
to zero and reverses the sign in the outer regions. Recalling that the worst perturbation 
is one of the orts of the subspace spanned by the columns of F, it can be expressed as a 
certain linear combination of the model functions. In this case, apparently, the differential 
scale terms [x 0]""" and [0 y]""" in combination with the radial distortion term [xr^ yi"^]"'" are 
the dominating contributors to the worst perturbation Un. 

How does one estimate the impact of the worst systematic error? Clearly, the term 
"worse" or "bad" can only be meaningful in a relative, or comparative sense. The singular 
values readily provide a meaningful comparison of various possible systematic errors. The 
most benign systematic perturbation is ui, and the ratio of the smallest possible to the largest 
possible systematic error is simply ak/cTi, which is the condition number of F. A very 
small condition number indicates an ill-conditioned problem. An inadeptly chosen distortion 
model may have a deficient rank, in which case one or more singular values equal zero, and 
the solution error can be infinite. For the JMAPS model in Eq. [21 the condition number 
is 0.14652, indicating a moderately conditioned model. The worst unit-norm perturbation 
results in a systematic error whose norm is 7 times that of the most benign unit-norm 
perturbation from the subspace of F. 



5. Orthogonal distortion models on the unit disk 



The two-dimensional Zernike functions (also known as the circle polynomials of Zernike) 
provide a ready-to-use orthogonal basis of scalar func tions on the unit di sk. These functions 
are used to decompose the optical aberration function (IBorn fc Wolflll999l ). and find other nu- 
merous applications in optics and engineering. Table [1] lists several low-order terms starting 
with Zq, which is a constant. An e xpanded set of the Zernike functions and transformations 
of algebraic monomials is given by iMatharl (120091 ). The Zernike functions are normalized to 
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-^/tt with a weight r: 

/ r ZlZlrdrdct> = 7i5mrAu (10) 
Jo Jo 

and each of the algebraic monomials in Eq. |2]is uniquely represented by a linear combination 
of Zernike functions. The polynomial model is not orthogonal because some of the terms 
include the same Zernike function, as can be seen in Table [H For example, the product 
of the 6th term (with ag) and 0th term (with qq) is equal to 7r/6. One can replace the 
non-orthogonal terms in model ([2]) with their significant Zernike counterparts, arriving at 

"Zr^ Z° ^Zf ^Z,-^ Zl 

_ ^z,-^ ^zf ^zf z-^ zr zl z\ 

This is the closest to the original model, which is orthogonal on the unit disk. It can 
be easily expanded to higher orders by including more terms constructed from the Zernike 
functions. 







a = 







(11) 



6. Orthogonal distortion models on the unit square 

Instruments of optical astrometry are equipped today with CCD detectors or other 
electronic sensors that are rectangular in shape. It is practical to use a functional basis, 
which is orthonormal on a quadrangle. Without a loss of generalization, a basis on the unit 
square can be applied to any rectangular field of view by means of coordinate normalization. 
The Zernike functions described in the previous section are not relevant, because some of 
them are not mutually orthogonal on the unit square. 

We demonstrate herewith, how a functional basis can be constructed on the unit square 
starting with a model, which includes polynomials or Zernike terms. This is implemented 
by the Gram-Schmidt orthogonalization process in the space of 2D scalar functions. If 
Ym are the original (non-orthogonal) vector functions, and Vm are the new basis functions 
orthogonal on the square, the algorithm is: 

Vi = Yi 

V2 = Y2-(Vi,Y2)/(Vi,Vi)Vi (12) 

V3 = Y3-(Vi,Y3)/(Vi,Vi)Vi-(V2,Y3)/(V2,V2)V2 

The inner products in this case are non- weighted integrals over the unit quadrangle, i.e., 

(V^,V„) = V^-V„dxdy. (13) 
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Normalized functions Vm = Kn/ a/ (Kn, Kn) have the additional advantage of all singular 
values being equal to unity for a dense and uniform set of sample points, which makes a 
perfectly conditioned model. 



A set of orthogonal functions Vj- 



m 



12, derived from the original JMAPS 



model ([T]) is given in Table [21 along with their norms on the unit square. The updated 
model of JMAPS distortions in terms of orthogonal functions is 



Vi V2 



V, 



V 



(14) 



As long as the field of view is well sampled with observations of calibration stars, the design 
matrix in these equations is nearly orthogonal, i.e., F"^F ~ 1^. 



7. Discussion 

Using the normalized vector fields based on the 2D Zernike polynomials for circular 
fields of view, Eq. [TTl or 2D orthonormal polynomials for rectangular fields of view, Eq. [TH 
leads to the most stable and well-conditioned solutions for instrument calibration parameters. 
Any systematic error of unit norm will result in a solution error of the same magnitude. A 
random error of unit weight (e = I, per Eq. E]) will result in an uncorrelated set of 
model coefficients, each having the same error expectancy, because the covariance of the LS 
solution is close to unity. This property of orthonormal distortion models comes in handy 
when we have to assess the significance of different terms in our model. On the one hand, the 
insignificant terms, which describe certain types of field dependent corrections that are not 
really present in the data, should be identified and eliminated as soon as possible, because 
the redundant degrees of freedom degrade the overall astrometric solution. This is especially 
important when the instrument is not very stable, resulting in numerous field-dependent 
calibration unknowns in the global adjustment. On the other hand, if the observational 
residuals prove too large, and their distribution suggests the presence of higher order terms, 
which are not captured by the distortion model, the orthonormal basis should be built up 
using the generalized Gram-Schmidt algorithm (^. 

Orthogonal plate models are also useful when possible internal degeneracies in the in- 
strument characterization have to be identified. A specific, but very common case is when the 
focal plane assembly includes more than one photodetector array, each being characterized 
by its own set of calibration parameters. In such hierarchical calibration models, lower-level 
parameters may be strongly correlated or completely degenerate with the higher-level pa- 
rameters. For example, a shift of each CCD in the same direction by the same amount is 
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indistinguishable from a common shift of the entire focal plane. Such double accounting of 
the same effect should be avoided, because it can lead to a complete solution failure, or diver- 
gence of iterated solutions. With an orthonormal plate model, the degeneracies are simply 
quantified by projection of the additional terms onto the previously accumulated basis. If 
the projections turn out to be significantly nonzero, the new terms can be orthogonalized 
analytically, as described in this paper, or numerically by the principal component method. 
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Table 1. Algebraic polynomials and Zernike function 



Zernike to Polynomials conversions Polynomials to Zernike conversions 



zl 



zr 

z,' 
zt' 



2r sin 
2r cos (f) 
v^(2r2 - 1) 
y/Qr"^ sin 20 
-\/6r^ cos 2(f) 
2v^(3r^ - 2r) sin^ 
2y2(3r^ -2r)cos(/) 



1 

X 

y 

+ y'^ . . . 

xy 

y' 



^0 

1 7+1 

2^1 

1 vO I 1 vO 
573^2 + 2^2 

1 /7-2 

1 7+2 I 1 7O I 1 7O 

i76^2 + i75^2 + 4^0 
1 7-1 I 1 7+1 

W2^3 +3^1 

1 7+1 I 1 7+1 
W2^3 + 3^1 
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Tabic 2. Vector basis functions orthogonal on the unit square 



Function 



Square norm 



Vi = 

V2 = 

V3 = 

V4 = 

V5 = 

V6 = 

V7 = 

V8 = 

V9 = 

Vio = 



1 



1 

X 




y 
y 




X ^ 

X2 - 1 

xy 

xy 

xr^ - if X 

2 14 

yr -Tsy 

5^2 I y2 _ 14 n 

4 

-gxy 

4 

-gxy 

x^ + iy^ + -il 



(Vi,Vi) =4 



(V2,V, 



(V3,V3) 
(V4, V4) 



(V5,V5) = | 
(V6,V6) = | 



(V7,Vr) 

(V8,V8) 
(V9,V9) = 



(Vio,Vio) = 
(Vii,Vii) = 



1984 
4725 



77344 
164025 



77344 
164025 



V 



12 



35^3 _ 27„„2 _ 30 

62 62 155 

35^3 _ 27 2 _ 30 

62-y 62 y 155-^ 



(Vi2,Vi2) = ^ 



